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Abstract 

Transient homogeneous nucleation is studied in the hmit of large critical sizes. Starting from pure 
monomers, three eras of transient nucleation are characterized in the classic Becker-Dor ing kinetic 
equations with the Turnbull-Fisher discrete diffusivity. After an initial stage in which the number 
of monomers decreases, many clusters of small size are produced and a continuous size distribution 
is created. During the second era, nucleii are increasing steadily in size in such a way that their 
distribution appears as a wave front advancing towards the critical size for steady nucleation. The 
nucleation rate at critical size is negligible during this era. After the wave front reaches critical 
size, it ignites the creation of supercritical clusters at a rate that increases monotonically until its 
steady value is reached. Analytical formulas for the transient nucleation rate and the time lag 
are obtained that improve classical ones and compare very well with direct numerical solutions. 
In addition, we propose and solve numerically a modified Becker-Doring model having a discrete 
diffusivity proportional to the area of a spherical cluster with k monomers for small k (as in the 
Turnbull-Fisher case) and to the cluster radius for large k (as in the case of diffusive growth of 
clusters). 

Key words: Kinetics of first order phase transitions; homogeneous nucleation; Becker-Doring 
equations; singular perturbation; asymptotic theory. 
AMS subject classifications. 82C20; 82C26; 34E15. 
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I. INTRODUCTION 



In first order phase transitions, it is possible to achieve that the system is in a metastable 
phase for values of the control parameter at which another phase is stable. Then nucleii of 
the stable phase appear at random locations and they typically grow or decay depending 
on whether their size surpasses a critical size 1]. In homogeneous nucleation processes, 
the probability density of finding a cluster of the stable phase with k monomers does not 
depend on the position. Homogeneous nucleation occurs in many examples of first order 
phase transitions, such as condensation of hquid droplets from a supersaturated vapor, glass- 
io-crystal transformations , crystal nucleation in undercooled liquids , and in polymers 
^, growth of spherical aggregates beyond the critical micelle concentration (CMC) p, 

1^ 



precipitation processes^ 
the miscibility gap 



and the segregation by coarsening of binary alloys quenched into 
In condensed systems, the time needed for the nucleation rate of 
supercritical clusters to reach a stationary value is large, which facilitates their theoretical 

and experimental study 31. Authoritative reviews of nucleation in condensed systems are 

n ^ n 

due to Wu [9] for theory, analytical and numerical methods, and to Kelton [3], who is 
not always reliable in his assessment of theories, but included a wealth of very valuable 
experimental data. 

While nucleation is a random process, the cluster size distribution function satisfies de- 
terministic equations which are particular cases of coagulation-fragmentation equations Q| . 
Typically, we assume that a cluster can grow or decay by adding or shedding one monomer 
at a time. Then homogeneous nucleation is described by the Becker-Doring (BD) discrete 



kinetic equations Isl llll. Suppose that nucleation occurs in a lattice in which there are 

nn 

many more binding sites, M, than particles, N, [y, We shall consider the thermody- 
namic limit, N oo with fixed particle density per site, p = N/M. Let pk be the number 
of clusters with k particles or, in short, k clusters, and let pk = Pk/M he the density of k 
clusters. Note that the number densities per site, p and pk, are both dimensionless. Number 
densities per unit volume are obtained dividing p and pk by the molecular volume, v = V/M. 
The BD equations (BDE) are 0] 

Pk = jk-i - jk = -D- jk, k>2, (1) 
jk = dk I pipk - pk+i \ ■ (2) 
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In (j2I), Sk is the binding energy of a k cluster, required to separate it into its monomer 
components. For spherical aggregates, 

Ek = ({k - l)a - - 1) ) ksT. (3) 
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This formula holds for /c S> 1, but we shall use it for all k > 1. aksT is the monomer- 
monomer bonding energy which, in the case of precipitation of crystals from a solution 
or segregation by coarsening of binary alloys, may depend on the particle density p (yolume 
fraction) through some empirical formulas In Eq. (jHl), o" = 27s(47rf ^/3)^/(A;b7'), where 
7s and v = V/M are the interfacial free energy per unit area (surface tension) and the 
molecular yolume, respectiyely. Note that a and a are both dimensionless. The correction 
3akBT/2 in Q ensures that ei = 0, and it irnproyes the agreement between the nucleation 



rate obtained from the BDE and experiments 



by Penrose et al 



More precise atomic models were proposed 



The monomer density pi can be obtained from the conseryation identity 

oo 
k=l 

in which the total particle density p is constant. In ((H), pk = dpk/dt and D±Uk = ±[uk±i—Uk] 
are finite differences. The time t, the discrete diffusiyity dk and the flux are nondimen- 
sional. t and dk are related to the dimensional time t* and decay coefficient dl as follows 
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t = nt*, 4 = §. (5) 

Here the factor Q has units of frequency, it depends on the particular model we choose for 
dk, and is determined in Appendix El for the Turnbull-Fisher (TP) kinetics (which assumes 
that a monomer has to oyercome an actiyation energy barrier for its transfer across the 
interface of a cluster). The TP dk is 

ion p-QKRT) 

dk = k'''e^-^^l\ ^ = -^^^, • (6) 

Here D = Dqc'^/'^^^^ is the diffusion coefficient in the liquid, Q is the actiyation energy 
for diffusion, R = ksNA is the gas constant and v is the molecular yolume. In the classical 
theory, dk is proportional to the surface area of a cluster. In other models, dk is selected 
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so as to yield the known expression for the adiabatic growth of a nucleus of critical size 



12| . The discrete diffusivity of these later models is 



by diffusion [7] or by heat transfer 
proportional to the cluster radius, thereby to k^^^. 

The flux jfc in size space is the net rate of creation of a A; + 1 cluster from a k cluster, given 
by the mass action law. Notice that we have selected the kinetic coefficient for monomer 
aggregation to dk, the coefficient for decay of a (fc + 1) cluster, so that 



is the equilibrium size distribution solving = 0. This is the detailed balance assumption 
whose validity is discussed in Wu's review |9i]. 

Eqs. dH), (j21), ©, Q and © form a closed system of equations that we can solve for 
an appropriate initial condition. If initially only monomers are present, we have pi{0) = p, 
and pfc(O) = for k > 2. The stationary solutions of the BDE and the phenomenon of 
phase segregation [3] will be described below in Section |Tll In Section IIIH we shall solve 
numerically the BDE for parameter corresponding to the glass-crystal transition in disilicate 
glasses PI and explain the results by means of an asymptotic theory valid for large critical 
size. A more detailed version of the material in Sections |n] and IIIII of the present paper 
can be found in Ref. which, in addition, contains a model for the kinetic constant dk 
based on thermal diffusion, not on activation processes as the TP model. The asymptotic 
theory for the discrete BDE is essentially the same for both models of the constant dk- 
Most previous asymptotic theories correspond to the continuum limit of the BDE which is 
a parabolic partial differential equation known as the Zeldovich-Frenkel equation (ZFE) Q], 
and thus differ from and are less precise than ours jl^ . 

A somewhat puzzling point in nucleation theory is that the growth of a cluster is 
proportional to its surface area (therefore to fc^/^) if we assume that capture or emission of 
a monomer is an activated process as in the TF model. However, clusters are supposed to 
grow by diffusion (and thus at a rate proportional to k^^^) after nucleation has stopped and 
coarsening begins jl^. In Section HVl of the present paper, we present a new modified BD 
model such that the growth of a cluster is proportional to k'^^^ for small size k (as in the 
TF model) and to k^^^ (as for the diffusive growth of precipitates) if k is large. We discuss 
briefly the results obtained with these two models and their relation to other studies in the 
literature. 
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II. BDE AND THEIR STATIONARY SOLUTIONS 



A. Equilibrium size distribution 

The equilibrium distribution ((Tj) satisfies jk = and it can be written as 

Pk = Pie-'\ (8) 
where Qk is the activation energy given by 

3 2 

gk = (Tk - (k - 1) (p, cTfc = - 1), (A; > 1), ip = \n {e°' pi) . (9) 

Here ai = = gi. Assuming k ^ 1^ g^ achieves its global maximum g^ = ok^J"^ 12 + 
akc^^^ — 3cr/2 at the critical size 

\ 3 

cr 



»: = fc.y. (10) 

Rewriting the flux (j21) in the BDEs in terms of the activation energy, we obtain 

3k = 4 { (e~^+^'= -I) Pk- D+pk} . (11) 

Eq. (PJ) is a spatially discrete Smoluchowski equation with diffusion coefficient dk and drift 
velocity 

Vk = 4 (e-^+^'^- - 1) . (12) 

Notice that < for k < kc (subcritical clusters shrink) and f ^ > for /c > kc (supercritical 
clusters grow). 

For the equilibrium densities (jH]), the conservation identity Q becomes 



oo 



e> = ^A;(eVi)'' e-"^ = ^fce'^^-"^ (13) 

k=l k=l 

This series converges for e"pi = e*^ < 1 (v? < 0), and diverges for e"pi > 1 [ip > 0). At the 
critical micelle concentration (CMC), pi = e~" {ip = 0), we obtain the critical density above 
which equilibrium is no longer possible. 



e°p, = l + ^A;e-"'=. (14) 



k=2 



For p > pc, the BD kinetic equations predict phase segregation, i.e., indeflnite growth of ever 
larger clusters, and there remains a residual monomer concentration whose density pic" — 1 



as t oo 



13|,|l4|. 
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FIG. 1: (a) Comparison of s„(i) evaluated (at different times) from the numerical solution of the 
discrete equations to the asymptotic result l|in|) (solid line), (b) K{T) calculated from Eq. 
H31|) with K(0) = (solid line) is compared to the numerically obtained position of the wave front. 
Data correspond to disilicate glass at 820 K. All variables are written in dimensionless units. 

B. The controlling parameters 

The simplest nucleation problem consists of solving the BD equations (@]), (Q) and (jllj) . 
with dimensionless activation energy gk = ctr — {k — l)(p, discrete diffusivity dk (to be chosen 
later) and initial conditions 

pi(0)=p, p2(0) = p3(0) = ... = 0. (15) 

The only parameters left in this initial value problem are p and a. p controls the long-time 
behavior of the BDE: If p < Pc given by ()14|) . pk{t) approach their equilibrium values (jH)), 
with monomer density pi that solves Equation li p > p^, cluster sizes grow indefinitely 
whereas their density becomes small. 

Let us identify the controlling parameters p and a in a physical system undergoing ho- 
mogeneous nucleation. A good experimental example for which abundant data exist is the 
transformation of lithium disilicate glasses to crystals (devitrification) . In disilicate, the 
free energy per molecule of the crystal phase in the activation energy (0) is proportional to 
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the undercooling 



ASKr„-r) 



whe.e r„ . the ,neM„g te.pe.atu.e, AS, is the .nola. ent.opy of fu.on a.d iV, . 

Avogadro's number; see the parameter values in Table I The dimensionless density 
p ^ g¥'(o)-Q can be extracted from Eq. ()16p as explained in Section IIIII In energy units, the 
activation free energy is ksTg^ = 7s47ra^ — ksTLpk^ where a = [3f /(47r)]^/^/c^/^ is the radius 
of a spherical k cluster. Thus 

ksT (^g,-^ + ^^= j,{4nY/\3v f'e/' - ASf{T^ - T)A;/iV^. (17) 
Comparing (fTTj) with yields a = {32nv'^ / 3)^^'^ 'js/ {ksT) , and the critical size 

The other parameters in Table I will be used later to model the discrete diffusivity in the 
BDE. We observe that the critical size increases with temperature: kc = 18 at 703 K and 
kc = 34 at 820 K. For other materials, such as undercooled liquid metals, critical sizes can 
be rather large: liquid iron at maximum und.ercooling has kc = 494, whereas kc = 2253 for 
liquid rutenium at maximum undercooling 3]. 

C. Equivalent Becker-Doring system 

As they stand, the BDE are rather stiff and hard to solve numerically This motivates 
the following change of variable 

p, = pie~^'su = e-"e'='^-'^'=Sfc, (19) 

according to 0- Note that = 1 in equilibrium. Since (71 = 0, this equation implies 

si = 1, (20) 

for all t. For the initial condition (fT3j) . e'^^^^~°' = pi(0) = p, and the conservation identity 
(jH) becomes 

00 

gv(0) ^^v^Y^f. e^^-^'Hk, (21) 

k=2 



Parameter 


Symbol 


Value 


Melting temperature 


T 

11 1 


1300 K 


Entropy of fusion 


ASf 

J 


40 J mol-i 


Surface tension 


Is 


0.15 J/m2 


Preexponential diffusivity 


Do 


2 xlO^ m^ s-^ 


Activation energy for diffusion 


Q 


440 kJ/mol 


Molecular volume 


V 


lQ-28 ^3 


TF time scale (703K) 




1.226 hours 


Critical size (703K) 


kc 


18 


Undercooling (703K) 


'f 


4.087 


Dimensionless surface tension (703K) 


a 


10.74 


Dimensionless free energy barrier (703K) 


n — Zh'^/'^ _ 3£ _|_ ,75 
i/m — 2 <^ 2 ' r 


25.177 



TABLE I: Data for lithium disilicate glass 
in which we have used (II 9j) . In terms of the Sk, the flux can be written as 

e"jfc = 4 exp[(A; + - ak+i] {sk - Sk+i), (22) 
and the BDE (H)) and ^ become 

s'k + Uk{sk+i - Sk) = -kifSk + dk-i {sk-i - 2sk + Sk+i), (23) 

for k > 2. Here, 

Ufc = 4-i-4e'^-''+"'=. (24) 

The term Uk D^Sk in Eq. ()23|) represents discrete advection, with a drift velocity Uk = 
—Vk + {dk~i — dk) ~ —Vki which is essentially minus the drift velocity in the original BDE 
for k ^ 1. Thus, the advection in Eq. (j^Hj) climbs up the activation energy barrier, from 
small values of gk to large ones. 

In summary, the transformed nucleation initial-boundary value problem consists of the 
balance equations (jSHI), the particle conservation equation (j2II), the boundary condition (f^ . 
Si = 1, and initial conditions 5^(0) = for all k >2. Its solution gives ip{t) and Skit) for all 
A; > 2 and all t > 0. 
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D. Stationary solution 



The stationary solution of the BDE has a flux independent of cluster size, so that e^jk = 
4exp[(/c + l)v?-o-/,+i] (sk-Sk+i) =j, from which (s^+i - Sfc) = -j exp[o-fc+i - (A; + l)v9]/4, 
and therefore 

Sk = l-Jl^ J , (25) 

1=1 ' 

for k >2. Since Soo = 0, j can be obtained from this expression in terms of an inflnite series 

Y.Zi exp[cri+i - {1 + l)v5 - In di] ' ^^^"^ 



Substituting back this expression into (j2SI), we obtain 



, _ 1 _ Eili^ exp[cr;+i - (/ + l)ip - \ndi] . . 



YaLi exp[o-;+i - (/ + 1)^5 - In di 
Then, = pie'^'^Sk- 



III. ASYMPTOTIC THEORY OF NUCLEATION 

In this section, we shall interpret the numerical solutions shown in Figures and |21 by 
using singular perturbation methods. Starting from pure monomers, the numerical solution 
of the BDE show that there are three well differentiated stages or eras of transient nucleation. 
After an initial stage in which the number of monomers pi decreases, many clusters of small 
size are produced and a continuous size distribution is created. During the second era, 
nucleii are increasing steadily in size in such a way that their continuum size distribution 
appears as a wave front advancing towards the critical size for stationary nucleation; see 
Fig. The nucleation rate at critical size is negligible during this second era. After the 
wave front reaches critical size, it ignites the creation of supercritical clusters at a rate that 
increases monotonically until its steady value is reached; see Fig. |21 Our asymptotic theory 
of the BDE will be described using the TF discrete diffusivity (jHl) and compared to numerical 
solution of the BDE for the crystallization of disilicate glass at different undercoolings. 



A. Initial transient 



Initially, pi(0) = p and there are no multiparticle aggregates. There is an initial transient 
stage during which dimers, trimers, etc. form at the expense of the monomers. This initial 
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stage is characterized by the decay of the chemical driving force = a + In pi to a quasi- 
stationary value (p, given by Eq. (fTBj) in the case of disilicate glass, and the emergence of a 
continuum size distribution. Knowing this, vje choose the initial chemical driving force (p{0) 
so that the quasistationary value (p given by Eq. is attained at the end of the initial 
stage. 

In materials such as disilicate glass at the temperatures we consider, the critical size is 
relatively small. Then v'(O) ~ ip, and the initial stage is very short. As the critical size 
increases (as in the case of undercooled liquid metals), <f{0) may differ appreciably from (p, 
and the initial stage lasts longer. However, even in such cases, the duration of the initial 
stage compared to the duration of the overall transient to stationary nucleation is of order 

fee''/' < 1 H- 

B. Wave front advancing towards the cluster of critical size 

After the first era, clusters of increasing size are formed. For sufficiently small clusters, 
the continuum size distribution approaches the equilibrium distribution with cp = cp. This 
situation can be observed as an advancing wave front in the variable Sfc(t), satisfying ~ 1 
(equilibrium) behind the front and ~ ahead of the front. This second era is described 
by Equations (pT|l to ()24|1 with ^ = (p and = 0. The critical sizes (fT^ for disilicate glass 
are relatively small, between 10 and 50, but they are large for undercooled liquid metals, 

— 1/3 

generally between 100 and 1000. Hence we shall use kc as a small gauge parameter 

Our asymptotic analysis will be carried out in the limit e — > 0, and therefore = e^'^ — > oo. 
Then dk, Uk and au in Eqs. (jHl), and (jSH) are smooth functions of /c > 0: 




(29) 
(30) 



u{k) = d{k - 1) - d{k) exp[v9 - a{k + 1) + cr{k)]. 



1. Position of the wave front 



In the numerical solutions shown in Fig. ^a), the graphs of Sk vs. k at fixed time have 
clear inflection points at some /c, where Sk ~ 1/2. The inflection point is taken as the 
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position of the wave front. In the continuum model, the front position k = kf{t) is a smooth 
function which obeys kf = u{kf). If we scale kf = K/e^, this equation becomes 
dK 

It 



U{K), T = et, 



(31) 



U{K) = \im[e^u{e-'^K)] = 2K^/^ sinh ^ ^{R-^^'^ 



1) 



(32) 



in the limit as e ^ 0. Fig. ^b) compares the wave front position calculated by solving ()3ip 
with K{0) = to the numerical solution of (|^. Note that the solution of (jHH) presents 
a time shift with respect to the numerical solution of the discrete model. This time shift 
reflects the breakdown of the continuum limit as i^" — 0, due to discreteness, and also the 
transient in ip{t) before it settles to (f. If the solution of Eq. (p?T|) - (p?^ is forced to agree 
with the numerical K{T) when the latter is, say, 0.1, the comparison fares much better. 



2. Shape of the wave front 

The leading edge of the wave front is a layer centered at K{T) in which Sk decreases from 
1 to as /c increases through it. The continuum representation of Sk in this layer is 

S{X,T;e), X = €^/^(k 



Sk 

Inserting (j^Hj) into (j^ . and then using (jHTj) . we obtain 



K 
Is 



i + f '(A-)A-|| ^ D(K) 



D(K) = lim 



e2 = ir2/3^Qgh(^(K-i/3-l 



(33) 

(34) 

(35) 



in the limit as e ^ 0. Had we carried out the same analysis for the ZFE, we would have 
found D{K) ~ d{e~^K) e. This would have resulted in a wider wave front and a longer time 
to ignition than those described below. 

3. Flux and wave front width 



Besides determining the shape of the wave front near its location, Eq. (j34p yields the 
behavior of the flux (creation rate of clusters larger than k) jk near k = kf. If we substitute 
© and into (j22I), we obtain 



Jk 



-1/2 ^2/3 g3<^/(2e) 



exp 



G{K) G'{K)X G'{K) G"{K) 



p3/2 



X' 



dS , , 
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Here, G{K) = (p {^K"^/^ — is a scaled version of the activation energy (^j). 

Since jk is proportional to dS/dX, it is convenient to differentiate (fH^ with respect to 
X in order to obtain an equation for J = —dS/dX, 



dT ^ ' dX ^ ' dX^ 



(37) 



Notice that J is locally conserved, and the following integral conservation identity holds: 

.00 



]oo 
— CXD 



dX 



dX 



JdX. 



Eq. has Gaussian solutions satisfying 

1 



J(X,T) 



exp 



X^ 



AA{T) 



(38) 



(39) 



which yields 



S{X, T) = - erfc 



X 



(40) 



for the wave front profile . Inserting Eq. (j39|) in Eq. pTjl , we find the following equation 
for A{T) > 0: 



^ - 2U'{K)A = D{K). 

After insertion of (jSHl), the flux (jH^ becomes 

G{K) G'{K)X G\K) 



(41) 



3k_ 

Joo 
Joo 



2A(p 



exp 



2e3 



,3/2 



2 ^ 



Que 



exp 



^ 3^ 
"2e3 ^ 2e 



X2 ,(42) 



(43) 



Here K = K{T) and A = A(T) are found by solving the differential equations (jHTj) and 
dH} with initial conditions K{0) = and ^4(0) = 3e^/(2(^), respectively. Eq. ^ is 
the classical Zeldovich quasi-steady nucleation rate of supercritical clusters, and it can be 
directly obtained from the stationary flux in the limit as e — > 0. 



C. The nucleation rate of supercritical clusters 



Let us now study the transient creation rate, in which j = j^^ increases from to the 
steady Zeldovich value As we have just seen, our theory predicts that the wave front 
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profile is given by ^Uj ), where K{T) and A{T) are solutions of Eqs. and \41\} , re- 

spectively. The flux of clusters with sizes larger than k is then given by Eq. Setting 
k = kc = (critical size) and X = (1 — K(T)) / e^^"^ in this equation, we obtain the 
nucleation rate predicted by our theory, j{t). Its integral over time yields the number of 
supercritical clusters, Nc{t). We shall consider now a different and more explicit approxima- 
tion of these results: We linearize the wave front equation 1)311) about the critical size K = 1 
and insert the solution into (021) • 



1. Linearization of the wave front speed about the critical size 



Let us fix /c = /cc = e ^ (critical size) in the definition ()HHj) of X: 



X='-^^K. (44) 



We now set X = K in ()42|) and perform the limit as e — > with k fixed. The result is 

J ~ j^e-^^y^-^'^'^^/^ ^ Jooe-^'''/^ (45) 



provided we use the limiting stationary value {4A)^^ = —G"{l)/2. 

The transient turns on when k = {1 — K)/e^^'^ = 0(1). Since f/(l — e^^'^n) ~ e^/^(^K/3, 
the wave front equation (jSU yields 

dn (p , . 

-T^ = --f^, (46) 

as e 0. The solution of this equation is 

n = KM e-^^'(^-^«)/3 = KMe-^'-'^'^/'-'^\ (47) 
r-' = l^e. (48) 

We select km as the value at which the flux j reaches its inflection point. Then we may con- 
sider that the wave front has ignited the nucleation of supercritical clusters. Straightforward 
use of Eqs. dHJ and (j46|l shows that 
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FIG. 2: (a) Evolution of the dimensionless flux at critical size j{t), and (b) number of clusters 
surpassing critical size Nc{t) as a function of dimensionless time for disilicate glass at 703K, kc = 18. 
Solid lines correspond to numerical results, dashed lines to the approximation given by Eq. (|42() . 
dot-dashed lines to the linearization approximation (|5U() and dotted lines to linearizing the equations 
for K{T) and A{T) as in Appendix C of Ref. 



Moreover, Tm = is the time to ignition, at which the wave front K{T) reaches the value 



K 



,3/2 



Km- From (jSH), we obtain tM and from (jlSj) and (HZj), we obtain the flux 
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-it-tM)/T 



Cm 



+ 



\-<?l'^KM 



In 



6e3 



+ 



Si^-^/s sinh [f (ir- 1/3 _ 1)] K-\ 



dK 



(50) 



(51) 



Here too is the duration of the initial stage. Integrating j{t) over time, we find the number 
of supercritical clusters as a function of time. In the limit as t — » oo, this number is 
Ncit) ~ joo (t — O), where the time lag 9 is approximately given by = t^/ + r7 + r£'i(e*"/^), 
in which 7 = 0.577215... is Euler's constant and Ei{x) is an exponential integral, see 
the derivation in Appendix B of Ref. The time lag 6 can be directly compared to 

experimental values j3[. 
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2. Comparison between different approximations 



Fig. |21^a) compares j{t) calculated from the numerical solution of the BDE for devitrifi- 
cation of disilicate glass at 703 K, from and with too = 0, and from Eq. with 
X = (1 - ir(T))/e3/2. We find that the more precise expression, Eq. (j42p . captures better 
the width and location of the transition region between j = and j = joo, as compared with 
the simple approximation given by Eqs. (j^UI) and ()51|) . Both approximations present a small 
overshoot and yield a smaller time lag 6 than that obtained from the numerical solution 
of the BDE. The overshoot decreases as the critical size increases. Another approximation 
consists of lineari zing the equations for K{T) and A{T) about the critical size = 1 as 
suggested in Ref. |l5|. This latter approximation is the worst one. This is not surprising as 
such approximation provides the same result for both the discrete BDE and the continuum 
ZEE. 

For disilicate glass at a lower temperature of 703 K, the critical size is smaller and our 
approximations deviate more from the numerical solution of the BDE, as shown in Figure 
|2I^a). Integrating j{T) over time, we find the number of supercritical clusters as a function of 
time, Nc{t), which is depicted in Figs. I21^b). The numerical solution of the BDE with the TF 
diffusivity yields a time lag 6 = 2.6. This value is close to those provided by the linearization 
approximation, 6 = 2.2, and by Eq. (jl^ . 6 = 2.3. Thus these analytical approximations to 
the numerical solution are reasonably good even for a relatively small critical size. However, 
6 = 2.6 gives 3.2 hours according to Table I, whereas the experimentally measured time 
lag is about 50 hours, cf. Fig. 5 of Ref. 0]. This discrepancy is due to having used the TF 
discrete diffusivity, which yields an excessively small time unit, as shown in Table I. 

IV. MODIFIED MODEL OF NUCLEATION AND GROWTH 

It is somewhat paradoxical that the TF discrete diffusivity is proportional to k"^^^ (cluster 
area), whereas growth stages of a cluster after nucleation are due to diffusive accretion of 
monomers which yields a discrete diffusivity proportional to k^^^ (cluster radius). In this 
section, we propose a model that interpolates between these two mechanisms and produces 
a discrete diffusivity proportional to fc^^^ for small cluster size and proportional to k^^^ for 
large cluster size. The idea is to consider that the discrete diffusivity should be consistent 
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with adiabatic growth of a large cluster at whose surface the concentration of monomers is 
different from the monomer concentration at infinity. We then modify the Becker-Doring 
equations to accomodate the resulting law for cluster growth Q]- 



A. Modified Becker-Doring model 

To be precise, consider a spherical fc-cluster which is growing adiabatically by diffusion 
of monomers across its interface. If the diffusion coefficient is approximately constant, the 
concentration solves Laplace's equation in spherical coordinates with densities pa at the 
sphere radius r = a, and poo at infinity: 

p(r) = Poo + ^° a. (52) 

Suppose now that the crystal particles occupy all available sites inside the cluster, whereas 
they occupy the volume fraction p outside the cluster. Then the number of molecules falling 
into the cluster surface per unit time is (1 — p)4:TTa^ {da/ dt) / v , and this number should equal 
the diffusive flux times the cluster surface area: 

Atto^D dp , . , 47ra^ da 

/ a = 1 - p — . 53 

V or V dt* 

This implies 

dk (poo - Pa) ^TraD {p^ - p„) {A7rf/^D{3ky/^ 



dt* il-p)v (l-p)f2/3 

Using now ()A9jl for large, near-critical cluster sizes, 

dk UDk^/'D+gk ,.^,,,,..2/3. ,-1/3 



(54) 



dt* t;2/3 



12D{k/vy^'{ak-'^'-^a), (55) 



in which ipa is the free energy per molecule at the cluster radius. Similarly, poo — Pa ~ '-p — ^a 
if the free energy per molecule is small, and Equations and ()55|1 yield 

,-1/3 (7r/6)^/^(^-^J ^(1 + ^) 

^"-^^ ~ (l-p)A:V3 ^'P'^^JV^r (56) 

Here we have used (j2Hl) and the definition 

(57) 
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Equations (PHj). (f^ and (fS7j) provide 

dt* l + efcV3/^ ' ^ ^ 

in which Q is as defined in (Q. 

Equation (j58|l describes the diffusive growth of a crystal nucleus in a glass phase such 
that the concentration at the interface is different from that at infinity. We shall now modify 
the BDE in such a way that the same equation is obtained in the limit of large clusters near 
critical conditions. We shall replace the BD fiux ((21) by: 

Jfc = (e-^+^'=+*'^Pfc - Pfc+i) 4, (59) 

which is equivalent to replacing pie^'' instead of the monomer density in We shall now 
select so that the new drift velocity 

Vk = (6-^+5'=+*'= - 1) 4, (60) 

becomes 

I 2/3 D+9k /„.,x 

in the limit of large, near critical cluster sizes. Writing time in dimensional units, Eq. ()6ip 
is dHHI) up to higher order terms. From © and (pljl . we obtain Vk ~ —k'^^^ — D^gk), 
which compared to (jU^ yields 

The modified BDE are therefore ^ with the fiux ()59p and (jU^ together with the mass 
constraint (Q. 

B. Equilibrium and stationary solution 

Our modified BDE do not have a Gibbsian equilibrium (jS|) nor satisfy detailed balance. 
Instead, the equilibrium distribution for which jk = is 

Pk = exp ( ~gk + 5^ j • (63) 
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As before, it is convenient to define a new distribution Sk which is 1 at equihbrium: 

Pk = exp ^-a + kip - ak + Y^ ^k, (64) 

for /c = 2, 3, . . . and Si = 1. In terms of Sk, the flux becomes 

jk = -4 exp ^-a + {k + l)ip - (Tfc+i + ^ j (sfc+i - Sfe). (65) 

The stationary solution satisfies the equation = e~"joo5 with constant flux j^o- The result 
is 

fc-i / I \ 

Sfc = 1 - joo X] I ^'+1 - + 1)V^ - *i - In rf/ j . (66) 
/=i \ j=i J 

Since — > as /c — > oo, we obtain the stationary flux 

joc = 1 ■ (67) 

ESi exp[ai+i -{l + l)ip- ^.^^ - Inrf,] 

C. Numerical results 



In terms of the Sk, the BDE are 



s'k + Uk{sk+i - Sk) = - [kip + ^^^1 j \ Sk + dk-i {sk-i - 2sk + Sk+ 



1), 



(68) 



for k > 2 and si = 1. Here the new advection velocity is 



Uk 



4-1-46'^-''+'^'=+*'=. (69) 



We have to solve the BDE with initial condition 5^(0) = for k > 2. The mass constraint 
(PTj) is now 



oo 



+ ^ A;e'=^-'^^+^i=i*^Sfc. (70) 



k=2 



In (l62|l . the parameter fi can be written as 



_ e(vr/6)^/^ 



because p = exp [99(0) — a]. 
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FIG. 3: Modified free energy = gu — Sj=i — a as a function of the size k for fi = 10 (solid 
line), = 1 (dot-dashed hne) and fi = 1/3 (dashed line). 

As indicated previously, with the initial condition of pure monomers (p rapidly evolves to 
a quasistationary value (f and a continuum size distribution emerges. We have to select the 
initial chemical driving force <f{0) so that the quasistationary value ip given by Eq. (fTI)|) is 
attained at the end of the initial stage. Notice that v^(0) enters the definition ()71|) and the 
constraint (f7n|) . Not knowing the monomer- monomer bonding energy a, we find easier to 
fix /i at different values and find if{0) in (f7n|) so that the quasistationary value (f of Eq. ()16|) 
is attained at the end of the initial stage. The limit as /i — oo gives us back the original 
BDE with the TP discrete diffusivity. 

We have solved the modified BD model with parameter values corresponding to disilicate 
glass at T = 703K {k^ = 18) in Table I for three values fi = 10, 1, 1/3. The corresponding 
modified free energy Ck = Qk — Yl^=i ^ j ~ « is depicted in Fig. El We observe that the 
energy barrier becomes flatter as /i decreases and the departure from the usual BDE with 
TF discrete diffusivity is greater. Fig. El^a) shows the evolution of the size distribution 
function Sk{t) towards its stationary profile when fi = 1. Comparison between the three 
stationary profiles of the size distribution function is depicted in Fig. |3I^b). Notice that the 
profiles for /x = 1 and /i = 1/3 almost coincide and are noticeably less steep than the profile 
for /i = 10. 

The size distribution function Sk has a step-like profile whose inflection point marks the 
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FIG. 4: (a) Profiles of Sk for = 1 at different times. The profile of the stationary size distribution 
function is indicated by circles, (b) Stationary profiles of the size distribution function for = 10 
(squares), = 1 (circles) and /u = 1/3 (dots). 




FIG. 5: Evolution of the wave front position towards the critical size kc = 18 for /x = 10 (dot-dashed 
line), fJ. = 1 (solid line) and ^ = 1/3 (dashed line). 

instantaneous location of the wave front advancing in size space. Fig. El shows the evolution 
of the wave front location. We observe that the velocity of the wave front decreases as n 
decreases. 

Lastly, we have calculated the evolution of the flux at critical size and the number of 
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FIG. 6: (a) Evolution of the normalized flux at critical size, j(t)/jooi and (b) number of clusters 
surpassing critical size, Nc{t), for disilicate glass at 703 K and ^ = 10 (dot-dashed line), ^ = 1 
(solid line) and fi = 1/3 (dashed line). Time lags are 2, 2.2 and 2.3, respectively. 

supercritical clusters for fi = 10, 1, 1/3; see Fig. IHl For these values of /x, we have found 
(^(0) = 4.79, 6.9, 21.8 and joo = 2.1 x 10-^ 4.1 x 10-^ 2.8 x lO'^, respectively. 

V. DISCUSION 

In this paper, we have studied the case of phase segregation resulting when p > pc. 
Previously, other authors had carried out asymptotic studies of the BDE in the simpler case 
of subcritical density, p < pc, in which initial conditions of only monomers, or more general 
ones, evolve towards the equilibrium distribution. In many cases of polynomial growth for 
dk, equilibrium is reached via a wave front profile for Sk, which is similar to Eq. (j40|) with 
A oc K^, and K oc T^, for appropriate positive 6 and p; see Ref. |l7| and references cited 
therein. This advancing and widening wave front leaves in its wake the equilibrium size 
distribution. 

In the more complex present case of phase segregation and indefinite aggregate growth, 
a quasicontinuum wave front of emerges after a short transient which is governed by the 
discrete BDE. After this, the leading edge of the wave front advances towards the critical 
size, and it slows down and stops there, leaving behind it a quasi-equilibrium state. The 
arrival of the wave front to the critical size marks the ignition of nucleation of supercritical 
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clusters, which ends when t 
theories (see 
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le stationary Zeldovich rate is reached. Previous asymptotic 
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22j | and references cited therein) have been derived for the 



continuum ZFE, not the discrete BDE, and thus their results systematically misrepresent two 
things: (i) the time lags for transient nucleation, as explained by Wu [9j, and (ii) the width 
of the wave front and the time to ignition in the nucleation rate. The latter discrepancies 
occur because the diffusion coefficient appearing in the continuum equation for the wave 
front satisfies Dbde{K) = Dzfe{K) — U{K)/2, and therefore the time to ignition in the 
nucleation rate for the BDE is smaller than the corresponding one for the ZFE. 

Let us briefl-V mention several existing asymptotic theories for the ZFE. Shneidman 



and Shi et al 



20| Laplace transformed the continuum ZFE and matched a first stage of 



pure advection of clusters to a local expansion about the wave front when it is near its 
final position at the critical size. They obtained our simplest formula for the nucleation 
rate, Eq. (fSUj) with the same relaxation time, ttf or ttdg, except that their values for tM 
were different from ()5H1 . This can be expected from Wu's arguments about approximating 



the discrete BDE by the continuum ZFE l^l; see the systematic shift of approximations of 
the ZFE with respect to numerical solutions of the BDE in Fig. 20 of Ref. l9|. Trinkaus 
studied a ZFE with a drift term linearized about the critical size (parabolic 
barrier) as an approximation to the full ZFE. Their results are comparable to those found 
yy means of the Laplace transform and matched asymptotic expansions; see Wu's review 
2|. All these authors obtained a transition region for the nucleation rate j{t) that was wider 
than observed in the numerical solution of the BDE. Several authors also found a nucleation 
rate for supercritical clusters that did not tend to joo as t ^ oo if A; 7^ kc j3, ll^ 
which is often called the asymptotics catastrophe 0]. Our theory is free from this deficiency 



12 1 . Shneidman 



22j | criticized Maksimov et al's result and extended his earlier asymptotic 



formula for the nucleation rate 



j to non-critical sizes. The previous criticism of using 
approximations to the ZFE instead of approximations to the discrete BDE apply to these 
works. Our more precise approximation using Eq. ()42|1 plus the exact equations for the wave 
front location and its instantaneous width improve upon other approximations and perform 
better for materials with large critical sizes. 

The time lag obtained from the numerical solution of the BDE with the TF diffusivity (or 
from our asymptotic approximations using it) is too small as compared with experimental 
results (about fifteen times smaller for disilicate at 703 K). The TF discrete diffusivity yields 
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an excessively small time unit, as shown in Table I. Another difficulty is that the TF dk is 
proportional to the surface area of the k cluster whereas diffusive growth of a cluster after the 
nucleation has ceased produces oc k^^^, as in late stage coarsening theories We have 
proposed a new modification of the BDE that yields a size advection velocity proportional 
to k'^^^D^Qk for small /c, and to k^^^D^g^ for large k, thereby interpolating between TF and 
diffusive growth formulas. With the new modification, time lags are somewhat larger, but 
joo increases unrealistically because the effective activation energy barrier diminishes. The 
usual way to improve agreement with experiments is to change the surface tension so that 
joo agrees with measurements. It is clear that we should optimize the choice of the surface 
tension and of the monomer- monomer bonding energy a to assess the merits of our proposal. 
Lastly, the flatness of the activation energy barrier for small /i as depicted in Fig. IHl indicates 
that expansions about the barrier maximum will provide poor approximations and that the 
corresponding asymptotic theory for our modified BDE may not be an immediate extension 
of our results for the classical BDE. 
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APPENDIX A: TURNBULL-FISHER DISCRETE DIFFUSIVITY 

To derive d^ for a glass-crystal transformation, we assume that the particles at the glass- 
crystal interface are undergoing reactions 

glass crystal. 

The forward and backward reactions have activation energies and During the con- 
version of one glass particle to crystal, one has to climb up one side of an energy barrier 
with height e+, and then descend an energy e_. The net free energy change should equal 
the chemical potential: 

= D^9k. (Al) 
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The forward and backward rate constants of the reactions are 

= (3e'T^, K- = /3e"^, (A2) 
respectively. Here /5 is a positive constant. The kinetic rate constants can be rewritten as 

K+ = /3e-^e-^+^*/^ K_ = /?e-^e^+^'=/^ (A3) 

provided Q/Na = {£+ + ^-)/2 is the mean barrier height. The common prefactor jSe"'^/^'^ 
is the inverse time constant of the reactions when Dj^gk = and neither phase is formed. 
For a spherical cluster with k monomers, this prefactor is the molecular jump rate at the 
cluster surface, which may coincide with the rate of molecular diffusion in the glass phase. 
Since a glass particle must push aside other glass particles when it moves, it is plausible 
that molecular diffusion in the glass phase is an activated process, as suggested by ()A3|) . 
The diffusion coefficient D is one third the mean velocity times the mean free path A ~ v^^^ 
according to elementary kinetic theory. Then the diffusion jump rate is SD/X^, and we can 
identify 

/S--*-^. fl = floe-*. (A4) 

Consider now the growth of a fc-cluster by converting glass particles to the crystal phase 
in the active sites at its interface. Individual conversions are discrete events whose time 
sequence is a Poisson process. Let fc+ be the number of glass to crystal conversions happening 
on the cluster interface in a short time 6t*. The expected value of is the number of active 
sites at the cluster interface times the number of glass to crystal conversions n^6t*. The 
number of active sites is the area of the spherical /c-cluster of radius [3kv / (An)]^^^ divided 
by the area a single molecule would occupy on the cluster interface, namely, 7r[3f /(47r)]^/^. 
Then, Ak = Ak"^^^ and satisfies 

(A;+) = AkK+5t\ {{k+ - {k+)f) = (A;+) = A^n+St*. (A5) 

The conversions from crystal to glass are another Poisson process for which 

{k_) = AkH-St*, {{k_ - {k_)f) = {k_) = AkH-St*. (A6) 

We assume that the backward and forward Poisson processes are independent, so 

{{k^-{k^)){k.-{k.))) = 0. (A7) 
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Equations ()A5|) to ()A7|) imply that the net number of glass-crystal conversions, 5k = k^ — k^, 
has the statistics 

{6k) = Ak{K+ - K^) 6t* , {{6k- {6k)f) = Ak{K+ + K^)6t*. (A8) 
Inserting ()A3|) and ()A4|) into these equations, we obtain 

- ^^'"i^ (^) m - mr-} - ^co=h [!^) St: (A9, 

We now identify the velocity Vk = {e^^+^'' — 1) dk in the BDE with Vl~^6{k) / 6t* , and this 
yields the diffusivity Eq. ()A9|) also provides the diffusion coefficient (j^ . 
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